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We present a new class of high-order imaginary time propagators for path-integral Monte Carlo 
simulations by subtracting lower order propagators. By requiring all terms of the extrapolated prop- 
agator be sampled uniformly, the subtraction only affects the potential part of the path integral. 
The negligible violation of positivity of the resulting path integral at small time steps has no dis- 
cernable affect on the accuracy of our method. Thus in principle arbitrarily high order algorithms 
can be devised for path-integral Monte Carlo simulations. We verify this claim is by showing that 
fourth, sixth, and eighth order convergence can indeed be achieved in solving for the ground state 
of strongly interacting quantum many-body systems such as bulk liquid 4 He. 

I. INTRODUCTION 

Many quantum Monte Carlo (QMC) techniques, such as path integral [ground state] Monte Carlo (PI[GS]MC) 
and diffusion Monte Carlo (DMC), rely on stochastic propagation of the Schrodigner equation in imaginary time. 
In all these methods, the probability distribution sampled is the matrix element, or the trace, of the imaginary 
time propagator 

G(r) =c- tH = c-^ t+v \ (1) 

with Hamiltonian H = T + V and kinetic and potential operators T — (—h 2 /2m) ^ V 2 and V — J2i<j v ( r ij)- 
Since G(t) is generally unknown, r is usually discretized into a sum of short time steps e so that the full prop- 
agator G(e) can be approximated by a product of short-time approximate propagator G(e). The computational 
effort, and therefore the efficiency of these QMC techniques depends on e. If G(e) is accurate to high orders 
in e, then a large e can be used to span a given imaginary time interval, resulting in fewer samplings of (but 
possibly more complex) G(e). 

For QMC simulations, there is a surprise lack of general higher order algorithms. For example, one has the 
well known second-order, primitive propagator 

G 2 (e) = e- eV/ / 2 e- £T e- ey / 2 = G(e) + 0(e 3 ). (2) 

(For computing the trace, any splitting first-order algorithm, such as e~ eV e~ eT , will also yield a second-order 
trace 1 -.) The highly successful pair density propagator, that approximates G by a pair-wise product of exact 
two-body propagators 2 - must also be second order in the general many-particle case, but possibly with a small 
error coefficient. However, this approach in practice is limited to solely spherically symmetric interactions due 
to the difficulty of evaluating the two-particle density matrix exactly. The only fourth-order method known for 
many years is the Takahashi-Imada 3 , Li-Broughton 4 propagator 

G TI (e) = e -^/2 e -eV-^ m[ V,[T,V]] e -eT/2 = G(g) + Q ^ > (g) 

where [V, [T, V]] — ^ J2i ^ iV\ 2 ■ This "corrector" propagator is only second order, but yields a fourth-order 
trace, as explained in RefjjJ. Thus until recently, there were only two second-order and one fourth-order algorithm 
for PIMC simulations. 

The problem of constructing higher order PIMC algorithms is the time-irreversible nature of the imaginary 
time Schrodinger equation. The short-time propagator can in general be approximated to any order by a product 
decomposition, 

N 
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with coefficients {U,Vi} determined by the required order of accuracy. However, in QMC applications, since 
(R'\ e~* ieT | R) tx e~( R A 415 ** 6 ) is the diffusion kernel with D = h 2 /2m, the coefficient ti must be positive in 
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order for the kernel to be normalizable as a probability distribution. As first proved by Shcng 5 and Suzuki 6 , and 
later by Goldman-Kaper- and Chin£, beyond second order, any factorization of the form (j4]) must contain some 
negative coefficients in the set {ti,Vi}. Thus, despite myriad of higher-order propagators of the single product 
form (j4|) for solving the time-reversible, real-time Schrodinger equation, none can be applied in PIMC beyond 
second order. It is only in the last decade that bona fide fourth order, forward algorithms with all positive 
coefficients have been found^^ and applied to DMC and PIMC simulation o 12 ! 13 ' 14 . In order to circumvent the 
Sheng-Suzuki theorem, one must include the operator [V, [T, V]] in the factorization process. Unfortunately, it 
not possible to go beyond fourth-order by including more operators. It has been showni& that a forward sixth- 
order propagator would have required the operator [V, [T, [T, [T, V]]]], which is non-separable and impractical to 
implement. More recently, by fine-tuning a family of fourth-order forward algorithm with two free parameters^ 
such that the fourth-order error is zero, Sakkos, J. Casulleras and J. Boronat 1 ^, and later also one of usi£, have 
achieved sixth-order convergence in computing the energy of a number of quantum systems including liquid 
4 He. Despite this spectacular advance, it must be noted that the fine-tuning must be done, in principle, for each 
individual observable. The algorithm is therefore only "quasi-sixth-order" rather than uniformly sixth-order. 

In this paper, we will present the first QMC simulations using a bona fide sixth-order and eighth-order 
algorithm for imaginary time propagation. The algorithm is based on the multi-product expansion 1 - of the 
short time propagator, which is a new way of achieving higher order convergence circumventing the Sheng- 
Suzuki theorem. This is reviewed in Section III Al followed by a brief introduction to path integral ground state 
Monte Carlo (PIGSMC) in section QTBl 



II. THEORY 



A. Multi-product expansion of G 

Let G2(e) denote the second-order split propagator then for a given set of n whole numbers {h}, the 
multi-product expansion of Ref|l9| yields the following second-order propagator: 

n 

G 2n (e) = J^dG^ie/k) = G(e) + 0(e 2n+1 ) (5) 
»=i 

where the expansion coefficient has the closed form 



n 



kf - ' 



(6) 



For PI(GS)MC, it is convenient to choose the sequences {ki} — {1, 2}, {1, 2, 4} and {1,2,3,6} to produce the 
following fourth, sixth and eighth-order propagators: 

&,(e) = -jG 2 (e) + jG^(i) (7) 

Gt(e) i a(e) + ^ai(i) (8) 
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As we will see later, these sequences are chosen because they are the minimal "commensurate" sequences. 
Schmidt and Lee2°- have previously suggested the use of (0 in path integrals and did use it in computing the 
two-particle density matrix. However, they did not suggest that it can be used for doing path integral Monte 
Carlo simulations. 

Since G(e) > 0, only the error terms in eq. ([5|) can be negative. Thus for sufficiently small e, these extrapolated 
propagators, despite the explicit subtractions, must be positive. Only when e is sufficiently large, the error terms 
overwhelm G(e) in a significant fraction of configuration space. We will see below that such large e cannot be 
used anyway because the propagators become highly inaccurate. One might argue that the error terms can be so 
singular that despite the smallncss of e, it can overwhelm G(e) at some specific locations. However, this cannot 
happen, because by construction G2(e) is bounded everywhere (except in the case of the Coulomb potential, 
which is a well known problem 4 - in PIMC and which we exclude from the present consideration). The subtraction 
of two bounded functions cannot be singular. This point will be clear when we present the explicit construction 
of extrapolated propagators and numerical results in the following sections. 
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B. Path integral ground state Monte Carlo 

The above multi-product propagators can be applied to any general PIMC simulations. Here, we will im- 
plement it in the specific context of PIGSMC. PIGSMC samples the whole probability distribution function 
corresponding to a discretized imaginary time propagation from a trial wave function *S? T (R) to the (in principle) 
exact ground state \&o(-R)> where R denotes all degrees of freedom, e.g. for the translational coordinates of N 
particles, R = (ri, . . . , rjv). 

For any trial wave function with non-zero overlap with the exact ground state, the exact ground state 
wave function can be obtained by evolving in imaginary time 

*o(#)oc lim [ G(R,R',^ T (R')dR'. 

/2^oo J 2 

G(f ) is evaluated by factorizing it into a product of small time step propagators G(e), e = ^jy, which can be ap- 
proximate by one of the above-mentioned short time approximations. Therefore, the full probability distribution 
to be sampled is 

P(Ro, ■ • • , R2m) =-j^^t(Ro) G(R , Ri;e) .. . 

G(R2M-l,R2M',c) ^t(R2m), 

so that the expectation value (^ol^l^o) of a local operator A(R) is evalutated by sampling A at the central 
time step, A(Rm)- For the energy, we take advantage of [H, G] = to obtain the energy estimator in terms of 
the local energy of the trial wave function El (R) = H^t / *t 

Eq = J dR . . . dR 2 M El(Ro) P(Ro, • ■ • , R2m) ■ 
These multidimensional integrations can be carried out with the Metropolis method. 



C. Implementing multi-product expansions in PIGSMC 

To see how one can implement these multi-product propagators in PIGSMC, we will now give a detailed 
discussion of the fourth order case. Considering G 4 at time step size 2e: 



G 4 ( 2e ) = l e -^ e -^ e -eV e -eT e -eV/2 



-e^ v e~ 2fT e-' v . 



In evaluating the matrix element of Gn(2t), since the first term on the RHS has one more operator e~ eT , it 
would require one more intermediate state integration than the second term, resulting in two dissimilar terms 
difficult to sample uniformly. The key contribution of this work is to enforce uniformity by artificially splitting 
the single operator e~ 2eT in the second term into two: 



G 4 ( 2e ) = t e -*V/2 e -eT e -eV e -eT e -eV/2 
o 

_ I e - eV e- eT e- eT e- eV , 



(10) 



which gives in the coordinate resprescntation 

(l|G 4 (2e)|3) = J d2(l|c- £T |2)(2|e- eT |3) 



/ 



3 3 
d2G (12; e)G (23; e ) e -^/2-eV 2 -eV 3 /2 F ( 123j e ) 



(11) 
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where we have denoted 14 = V(Rk) and abbreviated Rk — > k. We have defined 



F(123,e) 



V1+V3 



■v 2 ) 



and the free propagator 



G (12;e) = (l|c~ eT |2) = {^De)- m /'h 



(12) 



(13) 



We observe that: (i) Without the factor F, (fTTj) is just accurate to second order, (ii) By including F, only the 
potential energy is extrapolated in order to convert G to fourth-order. (Hi) For sufficiently small e, F > 0. (iv) 
If the potential function is mostly convex (such as Lennard- Jones type potential near the potential minimum), 
then one has 



V{Ri) + V(R 



(14) 



Since 



G (12;e)G (23;e) 



c 2 



H1 + H3 s2 



(27rDe) 3Ar / 2 



-G (13;2e), 



for fixed i?i and -R3, R2 is normally distributed about (Ri + Rs)/2 with width oc y/e. If i?2 is such that it is 
between R\ and R3, then the convexity condition (fl"4|) would guarrantee p2|) to be positive for all e. This only 
fails when the width of the Gaussian distribution for R 2 exceeds \Ri— R^\/2, suggesting that the near-positivity 
of F can extend over a rather wide range of e, which is indeed observed. Metropolis sampling requires exact 
positivity of F, which we ensure by using max(0, F), i.e. rejecting moves where F < 0. We also collect statistics 
about these rejections, so ensure that their rate is low, and decreasing with e. 
The generalization to higher order is now clear. For the sixth-order, Eq.®, 

G 6 ( 4e ) = ^ e -^/2 e -eT e -eV e -eT e -eV e -eT e -eV e -eT e -eV/2 

45 

_ i e -tV e -<T e -eT e -2eV e -eT e -eT e -eV 



1 

45" 



e -2eV e -eT e -eT e -eT e -eT e -2eV 



(15) 



yielding the coordinate representation 

G 6 (12345; 4e) = G (12; e)G (23; e)G (34; e)G (45; e) 

6^ -eVi/2-eVi-eV 3 -eVi-eVs/2 
45 

l„-eVi-2eV3-eV 5 , J_„-2 £ Vi -2eV & 

9 + 45 
Similarly for the eighth-order ©, 

G 8 (1234567; 6e) = G (12; e)G (23; e) . . . G (67; e) 



(16) 



54 



e -eVi/2 e -eV 2e -eV 3fr eV ie -eV 5e -eVe e -eV 7 /2 



L35 

27 
-40 £ 

, A -3eVi/2 e -3eV 4 -3eV 7 /2 J_ 

15 840 " 

For commensurate sequences one can factor out all the free-propagators and restrict the extrapolation process 
only to the potential energy function. 



- e -eVi e -2eV3 e -2eV6 e -eV7 



e -3eV le -3eV 7 



(17) 



III. RESULTS 



We have implemented the PIGSMC algorithm using multi-level sampling as described in the review^. We 
compare our new extrapolated fourth, sixth, and eighth order propagators with the primitive (second-order) 
and the fourth-order forward propagators 1 ' 4A 



G 4 A(e) 



e 6 e 2 e 3 



?[V,[T,V]] e -§T e - 



(18) 
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To demonstrate that our multi-product propagators work for realistic, and strongly interacting quantum 
systems, we apply them to the case of bulk liquid 4 He. We calculate the ground state energy Eq at equilibrium 
density p = 0.02186A" 3 , by a PIGSMC simulation of 64 4 He atoms in a simulation box with periodic boundary 
conditions. The decay time is (3 = 0.25K -1 , and we use the potential by Aziz et al.— . In Fig. Q]we show Eq/N 
as function of e for various propagators. We fit the polynomial a + be n (lines) to Eo(e)/N, where n is the order 
of the respective propagator. Since the order of Eq{e) is defined as the e — > behavior, we have restricted the 
fits to small values of e - the end point of the lines indicate the fitting interval. These propagators are compared 
at equal time steps: G 2 (e), G±(c), G 6 (e), and G 8 (e), to verify the order of convergence. 

The primitve second-order propagator (open circle) is clearly a poor approximation, with a large error even 
for small e, and therefore requires a large number of beads. The simplest fourth-order forward propagator 
4A, Eq. (|18[) . is a significant improvement, as can be seen in the behavior of E (e)/N (open square), with error 
coefficient smaller than our fourth-order multi-product propagator (jlip (filled square). However, the forward 4A 
propagator requires the the computation of [V, [T,V]] cx |Vil^| 2 , and its relative efficiency would depend on the 
complexity of evaluating this gradient. Both can be fitted well by a fourth-order polynomial with n — 4. Finally, 
the closed triangle and circle show the convergence of the sixth order (fT6|) and eighth order (TiTl) multi-product 
expansion, which indeed has a smaller e dependence in the range of Fig. [1] These multi-product propagators are 
true high order propagators and will produce sixth and eighth order convergence for the expectation value of 
any observable. The present results constitute the first implementation of a quantum Monte Carlo simulation 
with a bona fide imaginary time propagator of higher than fourth order. At small values of the time step size, 
say at e = 0.005, the sixth and eighth order algorithms produce very precise results which are not indicated by 
Figffl 

The multi-product expansion of G, Eq. ([5]), is not strictly positive everywhere for a finite time step e, as we 
have discussed above for the fourth-order case. In Fig. [2] we show the ratio R n of attempted MC moves that 
are rejected due to negativity of the multiproduct expansion. R n is decreasing with e as expected. Only in the 
sixth-order case that we observed a non-monotonous behavior at very large e, that the ratio R n decreases with 
increasing e. This occurs at large values of e where the error of the energy Eo(e)/N is rapidly increasing with 
e (outside the range of Fig. [!}. This may due to system configurations with complicated intertwining negative 
and positive regions of G. We want to stress that, since this happens only for e much too large for quantitatively 
correct results, it poses no practical limitations. 

The convergence plot Fig. [1] does not reveal the actual computational effort required for a desired accuracy. 
From our derivation, it is clear that the computational effort of G4, Gq, and Gs are roughly equivalent to 
running G2 twice, four, and six times, respectively. This then means that for a given e for G2(e), one should 
compare it to G4 at 2e, G§ at 4e and Gs at 6e, that is, for an equal effort comparison, we should compare G2(e), 
G4(2e), Ge(4e) and Gs(6e). This is done in Fig. [3J In this comparison, at a given e each algorithm uses the 
same number of beads. For the forward algorithm 4A, this comparison neglects the additional cost of evaluating 
the gradient |ViT^| 2 . If |ViT^| 2 required no more effort than that of evaluating the potential, then propagator 
4A would actually outperform the extrapolated sixth- and eighth-order algorithms at time steps e > 0.003K -1 . 
This confirms that, in principle, a purely forward time step algorithm can be more efficient, provide that |ViV| 2 
can be easily evaluated. However, the higher order extrapolated algorithms are clearly easier to derive and 
implement. Morever, when very high accuracy is required, such the "chemical accuracy" required in quantum 
chemistry applications, then a higher order algorithm will always outperform a lower order algorithm. This is 
specially critical in determining the equilibrium configuration or conformation of clusters and macromolecules, 
where energy differences are very small. 

IV. CONCLUSION AND OUTLOOK 

In this work, we have shown how to implement the multi-product expansion of the imaginary time propagator 
in QMC for solving strongly interacting quantum many-body systems, such He, to any desired order in the 
imaginary time step e. In particular this work is the first demonstration of truly sixth- and eighth-order QMC 
algorithms. In the case of 4 He, our results suggest that these higher than fourth-order algorithms may not be 
more efficient than purely forward time step fourth-order algorithms, but they do have the simplicity of not 
requiring the potential gradient. This is particularly useful in simulating non-cartesian coordinate systems, such 
as molecules 2 ^ with anisotropic constituents and rotational degrees of freedom. Moreover, these extrapolated 
propagators are the only higher order algorithms possible in cases where the double-commutator cannot be eval- 
uated, such as for the diatoms-in-molecule potential 2 ^. Finally, for QMC applications where chemical accuracy 
is required, such as in determining equilibrium configurations and conformations, our sixth and higher order 
multi-product propagators will be computationally more efficient than fourth-order propagators. 
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FIG. 1: (color online) Ground state energy Eo of bulk 4 He, simulated by 64 4 He atoms, as a function of imaginary time 
step e. Decay time was j3 = 0.25K -1 . We compare results produced by the primitive second-order propagator G2(e) and 
the fourth-order forward propagatopi^ G4a(c) (denoted "4A") with our fourth, sixth, and eighth-order multi-product 
propagators Gi(e), G§{t) and Gg(e), (denoted "MP"). Each Eo(e) is fitted with the appropriate polynomial. 
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FIG. 2: (color online) The ratio R„ of rejected MC moves that would lead to a negative propagator G. R„ is decreasing 
with time step e. 




FIG. 3: (color online) A roughly equal effort comparison of algorithms (?2(e), Gi{2e), Ge(4e), Gs(6e) and G4a(2e) for 
the same ground state energy Eq as in Fig. [T] 



